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Abstract 

We present a numerical method to deal efficiently with large numbers of par- 
ticles in incompressible fluids. The interactions between particles and fluid are 
taken into account by a physically motivated ansatz based on locally defined drag 
forces. We demonstrate the validity of our approach by performing numerical simu- 
lations of sedimenting non-Brownian spheres in two spatial dimensions and compare 
our results with experiments. Our method reproduces essential aspects of the ex- 
perimental findings, in particular the strong anisotropy of the hydrodynamic bulk 
self-diffusivities. 
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I Introduction 



An important paradigm problem of multiphase flow, incorporating aspects of both a 
discrete particulate and a continuous liquid phase, is the sedimentation problem. Here 
the particulate phase slowly settles into a packed bed at the bottom of a container filled 
with viscous liquid. The main theoretical challenge in the sedimentation problem is the 
understanding of the long range liquid mediated hydrodynamic interactions between the 
immersed particles. Many aspects are well studied and understood, as are for example 
the mean sedimentation velocity (y) in the Stokes flow regime for a monodisperse 
^ 1^ or a polydisperse system 0]. Less well understood are fluctuation phenomena, for 
example the properties of the particle velocity distribution, as, e.g., its variance, time 
autocorrelation, and effects arising from particle polydispersity. Recently, much attention 
has been focused on the long time behavior of the particle motion as described by strongly 
anisotropic hydrodynamic diffusion constants 0. Moreover, most analytical work pertains 
to the regime of dilute suspensions, whereas semi-dilute and concentrated suspensions still 
present a challenge to theoreticians. 

The focus of this paper is the application of a novel computer simulation method 
for particulate multiphase flow ||^ to the sedimentation problem. Among the standard 
computational approaches to multiphase flow in general and the sedimentation problem 
in particular are (i) continuum descriptions 0, ^ that neglect the discrete nature of the 
particulate phase, (ii) methods that exploit the linearity of the Navier-Stokes equations 
at very small Reynolds numbers for either the simulation of hydrodynamically interacting 



point Stokeslets [T0| or extended spheres |TT|, |T^, |T3[, (iii) finite element and finite 
difference techniques in two (2D) or three dimensions (3D) that represent the proper 
Navier-Stokes equations with no-slip boundary conditions for the liquid on the particle 
surface, and (iv) algorithms that combine a continuum discription for the liquid phase 
with a discrete representation of the particulate phase . 



All of these approaches have inherent strengths and weaknesses. The continuum ap- 
proaches (i) suffer from difficulties to represent the particulate phase in terms of continuum 
quantities as pressure and stress. Both the Stokes flow methods (ii) and the differenc- 
ing schemes (iii) give very reliable results where they apply. However, the Stokes flow 
methods (i) are in principle restricted to low Reynolds number flow and inherent to the 
differencing schemes (iii) is a tremendous computational effort that limits their appli- 
cability to systems with very few particles ||15|. We therefore choose a method of class 
(iv) comprising basically a Lagrangian integration of particle trajectories coupled by local 
drag-force terms to a Eulerian Navier-Stokes integration. 

Such a method allows immediate access to basically all physically relevant quantities 
in the system, including particle coordinates and both particle and liquid velocities, at 
little more computational costs than a standard real space Navier-Stokes integration. The 
main drawback is that a neglect of the proper boundary conditions between particles and 
liquid will result in inaccurate rendering of the short scale flow properties. However, we 
will be mainly concerned with collective phenomena, i.e., the effects that arise when the 
number of particles is large and we do not have the ambition to describe accurately the 
local flow fields [|l^ on the scale of the particle size. Single particle and particle pair 
movement can to a high degree of accuracy be treated analytically |T6[. 

Nevertheless we want to see which and how phenomena emerge directly from simple 
modeling assumptions — as opposed to using semi-empirical expressions as, e.g., done 



in Refs. [14, M. In particular, we rely on the fact that the long-range hydrodynamic 



interactions, which we presume to be the most important for collective phenomena, are 



9. 



correctly represented by the velocity and pressure fields evolved by the Navier-Stokes 
integration. 

In this paper, we will first describe a two-dimensional (2D) implementation of our 
simulation technique in Sec. ||. We will then describe simulations of the sedimentation 
problem in a broad range of particle solid fractions. We first show particle trajectories 
resulting in the sedimenting bulk (Sec. |III.A|) . As function of the solid fraction, we 



address and discuss the hydrodynamic diffusion (Sec. [lll.BD , the hindered settling function 
(Sec. IIII.C]) and particle velocity fluctuations (Sec. [1II.D|). 



II Simulation Technique 
II. A Modeling assumptions 

Our aim is to capture the large scale collective processes of the combined particle-liquid 
system with a computer simulation. Therefore we need to model comparatively large 
systems in moderate time. We aim at the development of a tool that allows simple checks 
of qualitative and quantitative behavior of multiphase flow. To this end we concentrate 
on the momentum exchange between particle and liquid phase. Our simulation technique 
makes physically motivated assumptions to simplify the computational difficulties asso- 
ciated with a full treatment of the no-slip boundary conditions. Our assumptions are as 
follows. 

(i) For small Reynolds numbers lubrication forces dominate as particles approach 
each other or the walls of the container. These forces of dissipative nature theoretically 
prevent particles from touching each other. Since we are interested in long-range effects, 
we neglect the short-range lubrication forces, and use instead a simple contact model 
between particles. 

In 2D, the particles i are represented by disks with radii chosen from a Gaussian 
distribution hp{r/f) with mean f, which is cut off at its standard deviation a, 



hp{r/r) oc 




1| <p. 



Here, the distribution is written in terms of the dimensionless radius r/f and the poly- 
dispersity parameter p, which we define for the purposes of this paper to be the standard 
deviation of the original, full Gaussian divided by the mean radius, i.e. p = a/f. 

Our contact model assumes elastic, central, and repulsive forces, which act when two 
particles i and j touch each other M, O, IT^, IT^ . That is, the forces between two particles 



are zero unless the distance between their two centers is less than the sum of the radii 
Ti + Tj. Then the acting force is assumed to be 

Here Xj = {xi,yi) is the position vector of the center of particle i and fiji is the unit 
vector pointing from the center of particle i to the center of j. The spring constant fc„ 
(cf. Table |l|) describes the repulsion due to lubrication forces and depends on the size 
and shape of the colliding bodies. Tests including a velocity proportional damping term 
in Eq. (^ show that our simulation results do not significantly depend on the details of 
the particle-particle interactions. 



Apart from the drag forces F^, which will be addressed below, we further add gravity 
(including buoyancy) as a contribution to the total force acting on each single particle, 

^ ^ 4 

i^tot = Fd + Fei + -Tirfipp - pi)g. (3) 

Here, pp and p\ are the densities of the particulate and the liquid phase respectively, and 
g the earth's gravitational acceleration which is assumed to act in —y direction for the 
rest of this paper. Once we know all forces, we employ a fifth order Gear predictor- 
corrector integration algorithm to obtain the particle trajectories [^. Rotational degrees 
of freedom are not considered. 

(ii) We will assume that the exact arrangement of particles and the exact distribu- 
tion of the interstitial liquid on short scales do not affect the qualitative features of the 
large-scale flow pattern. As already mentioned above this means we neglect detailed lu- 
brication forces resulting from the no-slip boundary conditions of the stress tensor at the 
liquid-particle interfaces, replacing them by the above described particle interactions. For 
particle Reynolds numbers Re = rpi\Vi — Vi\/rj smaller than one, we propose to model the 
momentum exchange between single particles and the surrounding liquid by means of a 
Stokes-type drag force, 

F^ = CdVn{V,-vl), (4) 

here vl denotes the velocity of particle i, is an average local liquid velocity at the 
position of the particle and t] denotes the liquid's kinematic shear viscosity. Equation 
with Cd = Qir holds exactly in three dimensions (3D) for a single sphere moving with 
respect to the liquid resting at infinity. We obtain the Vi by linear interpolation of the 
liquid velocities from the four grid points of the quadrilateral MAC mesh of the liquid 
integration surrounding the particle center. 

In our case, is a parameter. For all choices of grid and particle size we determine 
Cd such as to obtain the Stokes velocity for the simulation of a single sphere falling in 
a viscous liquid. Commonly^, we use grids sufficiently large so that does not differ 
from Qtt by more than a factor of 2. The main physical effect of this local drag force is to 
introduce a tendency to locally equalize the liquid and particle velocities, similar to the 
no-slip boundary condition that implies particle and liquid velocities being equal on the 
particle surface. 



II.B Fluid Model 



To describe the time evolution of the velocity and pressure fields v and p of the fiuid in 
our simulations we start with the Navier-Stokes equations. 



Pi 



dv 

at 



+ {vV)v 



(5) 



Here, pi denotes the liquid density and r] its shear viscosity, / is the volume force den- 
sity in the system, including gravity, drag-force density and buoyancy response. Since 
we have in mind applications to liquid systems or gaseous systems with typical veloc- 
ities much smaller than the velocity of sound, we complement these equations by the 
incompressibility constraint, 

Vv = 0. (6) 

^ Wc find that the "drag coefficient" Cd depends on the dimcnsionless ratio of particle radius to grid 
spacing (see Sec. II.B), r/ Ax. In the regime of particle Reynolds numbers Re w 0.5 for our simulations 
we have determined the empirical relationship C^/GTr = 1 + 6.27{Ax/r)~ 



-1.5 
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Equation (Bf) presents a constraint on the velocity field that must be fulfilled at all times 



and may be employed to obtain the pressure field via an iterative procedure |22|. We will 



sketch the basic ideas briefiy and refer the reader for more details to Refs. [pT|, |23] and 



23| 



The pressure as well as the velocity components and Vy are located on three quadri- 
lateral meshes with lattice spacing Ax. With respect to the pressure grid, the grids for 
the X and y velocity components are shifted by Ax/2 in x and y direction respectively 
— this construction is commonly referred to as the MAC mesh and has several computa- 
tional advantages, i.e., it is a simple means to avoid numerical instabilities due to mesh 



decoupling 



We obtain the pressure and the velocity components by an iterative procedure. Let 
the index n refer to values at time t = and n -|- 1 to those at t = t^+i = t„ + At 
after a timestep of duration At. The index k shall denote an iteration index. We define 
Pn+ifi = Pn, i-e., we start an iteration for the new pressure at time t„+i with the old values 
at tn- We obtain a tentative velocity field at t = tn+i from an evaluation of the discretized 
Navier-Stokes equation, 

^^ ^n+l,fc^l Vn ^ _ + /„ + r]V'^Vn. (7) 

A von Neumann stability analysis ||2l], |2^ of the described discretization of the Navier- 



Stokes equation shows that the values Ax and At are subject to the two stability con- 
straints 

4?7 



and 



At < ^ 9 

In general, the velocity field Vn+\,k^\ resulting in does not satisfy the incompress- 
ibility constraint (P). To the end of lowering the modulus of the divergence of the velocity 
field — in so far it differs from zero — we calculate a new, likewise tentative, pressure 
field 

Pn+l,fc+l = Vn^l.k - ApiVtT„+l,A;+l. (10) 

Here, a large value of the parameter A is crucial for rapid convergence. It is, however, by 
stability requirements constrained to 

0<A<i^. (11) 

- - 4At ^ ' 

Subsequently new velocity values consistent with the may be obtained using 

Eq. (|^, or equivalently, without reevaluation of the Navier-Stokes equation, using the 
relation 

Pi = V(p„+l,fc+l - P„+l,fc). (12) 

Iterating Eqs. (p!oD and ([T^) yields a new pressure field and after convergence a con- 
sistent, divergence free velocity field for time tn+i- This iteration is equivalent to solving 
a Poisson equation for the pressure. 

For our purposes, the described algorithm has three advantages. It (i) generalizes 
straightforwardly to 3D and (ii) unlike spectral or streamfunction methods it gives im- 
mediate access to the quantities p and v in real space. Fast access to the latter is crucial 



for the drag-force calculation in Eq. (^). Moreover (iii), only the chosen coarseness of 
the spatial and temporal discretization limits the range of Reynolds numbers addressable 
in the simulation. The use of multi-grid techniques instead of the above local iterative 
method will improve the speed of our method for large system sizes. 

For the present paper we have implemented the above set of equations in 2D. The 
uncoupled liquid and particle equations of motion are directly obtained by dropping the 
third 2;-component of the velocity and applying 2D versions of the differential operators 
to the liquid equations. The coupling of the two sets of equations, however, presents some 
problems. A 2D simulation should, strictly speaking, be one of rigid, parallel cylinders 
moving in a liquid that only displays motion in the xy plane. However, at a fixed small 
Reynolds number, the drag per unit length of a cylinder |2^ does not depend on its radius 
in an infinitely large system. Thus, qualitative effects of particle polydispersity cannot 
be expected to be present in such a model. We therefore chose to refer the Stokes drag 
Eq. to a reference length z equal to the average particle diameter. Alternatively, one 
may think of z as the depth of the liquid filled vessel, in which 3D spheres move such 
that their motion is constrained to lie in the xy plane. The Stokes drag per length then 
enters the particle and liquid equations as coupling term. Since we have not included 
any Brownian motion of the particles our model is only applicable in the non-Brownian 
regime of large Peclet numbers Pe ^ 1, characterized by comparatively large particles 
and low temperatures (cf. Sec. [1II.D| ). 

All our simulations are performed in a rectangular vessel of width and height Ly. 
On the walls of the container we impose no-slip boundary conditions, i.e., the liquid 
velocity there is taken to be zero. The particles are initially placed randomly and without 
overlap within the vessel. Their size is characterized by their mean radius f and the 
degree p of polydispersity, which is kept constant and equal to 0.04 in all the simulations 
presented in this paper. See table |l| for a detailed list of simulation parameters.^ 

In the simulation the particles are point like and their volume only comes into play 
via the drag force. Presently, the model does not take into account the fractions of solid 
and liquid in a given volume element. This fraction, which is very important to describe 
backfiow, will be included in a more extensive study p. 



Ill Results and Discussion 

In the remainder of this paper, we wish to show that the set of equations above — 
even if they might appear to be rather crude — constitute a reasonable description of a 
particulate two-phase system, and are, in particular, capable of reproducing the essential 
aspects of a sedimentation process. 



III. A Particle trajectories in batch settling 

We first consider the particle trajectories in a batch settling experiment. Typical ex- 
periments |25| examine the behavior of a mixture of glass beads in a liquid kept in 
a container. The glass beads start to sink leaving behind clear liquid at the top and 
particle layers on the bottom of the vessel. Index matching techniques allow to trace 
the trajectories of single particles in the mixture H]. In the limit of infinite dilution, in 



^ A typical simulation run is characterized by about 3 200 nodes for the liquid integration and about 
1 500 particles. It involves about 10'' full timesteps for the liquid and 10^ timesteps for the particle 
integration. On a SPARC 20 workstation the required CPU time is of the order of lOh. 



sufficiently large vessels these trajectories are straight lines from the top to the bottom of 
the container. The particle velocity is the Stokes velocity Vs{ri) = {2/9)grf{pp — p\)/rj, 
the terminal velocity of a sphere of radius settling in an extended viscous liquid under 
the influence of gravity. For the purposes of this paper we will neglect the influence of 
polydispersity and set Vs = Vs{f). 

At finite concentrations, particles start to feel each others presence and the trajectories 
are no longer straight. One observes very complicated collective effects of the suspension 
which organizes into continuously changing regions in which particles move in different 
directions. Denser regions with more particles per unit volume drag liquid with them and 
the displaced fluid travels upward through less dense regions that present less resistance 
to the flow. On smaller scales particle velocities change due to close encounters. The 
resulting typical particle trajectories are highly irregular. Even in the laboratory frame, 
one finds particles moving upward in the suspension, describing loops and displaying 
settling tendencies only if observed for long enough time. 

For illustration we show in Fig. |l]a the trajectories of 20 randomly selected particles 
out of a 1000 settling in a suspension at solid volume fraction $ = 0.157. Here, the solid 
fraction is the dimensionless ratio of the area in the xy-plane covered by disks to the total 
area in the simulation, i.e., 

^^J2^^'/L.Ly. (13) 

i 

In Fig. [T|b we display the same trajectories as viewed from the comoving frame, i.e., we 
subtract the mean settling velocity from the particle velocities in the lab frame and shift 
the trajectories such that they all start in one point. For long enough times each trajectory 
displays the characteristics of a random walk. Its direction changes frequently on scales 
larger than a coherence length, which measures the average distance that a particle has 
to travel before its motion becomes independent of the initial velocity. 

In fact, experiments find that the particle velocities at large times decorrelate 0. 
However, the motion is anisotropic, spreading much more rapidly in the main settling 
direction than perpendicular to it. 

III.B Diffusion constants 

One may describe the long time behavior of the particle displacements in a sufficiently 
large system as a diffusion process in the frame comoving with the settling particles. The 
diffusion is characterized by its mean square displacements in vertical y and horizontal x 
direction, 

(At)) = ^Ek^W-^^(o)f , (14) 

i 

% 

The sum runs over those particles that are above the final bed height for t ^ oo at the 
given moment t. From their initial positions we determine the average motion to obtain 
{V). 

As already indicated by Fig. |T]b, a prominent feature of hydrodynamic diffusion in 
batch sedimentation is the strong anisotropy of the process with different diffusion con- 
stants Dx and Dy in horizontal and in vertical direction. Both Dj. and Dy are obtained 
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from the long time behavior of the particle displacement by the relations 

{x\t)) ~ 2D^t, (16) 
{y\t)) ^ 2Dyt. (17) 

In practice, the squared particle displacements saturate due to the finite system size 
after passing a regime linear in time. In this regime, one can determine and Dy 
using Eqs. (|16| and |l^. We show typical example for the time dependence of the mean 
square displacements in Fig. |^. We have normalized both x and ^-displacements by 
division through the squared particle radius, i.e., {xD = (x^(t))/f^ and accordingly {yl) = 
(?/^(t))/f^. We obtain the dimensionless time t* by division of t through f/Vs- The 
straight lines superposed to the data shall emphasize the approximately linear behavior 
for intermediate t*. The slopes of these lines, averaged over several runs with different 
initial particle positions, are equal to twice the dimensionless diffusion constants D* and 
D* which result from division of and Dy by the quantity fVs- 

In Fig. ^ we plot D* and D* vs. the solid fraction $. As is observed in real 3D 
experiments, also here the diffusion constants increase for small $ until in the simulations 
a maximum at area fraction $ ^ 0.12 is reached. The location of the maximum in the 
experiments is at volume fraction $ f» 0.1 at a slightly lower value. At increasingly 
higher $ one observes a renewed decrease of the diffusion constants. The values of our 
2D diffusion constants lie by a factor of ~ 2.5 ... 3 above the 3D experimental values ||^ 
or 3D simulation results using Stokesian dynamics methods |]I2|, |T3|. 

Qualitatively, we understand the existence of a maximum at intermediate volume 
fractions. For $ 0, hydrodynamic interactions between the particles are not very 
pronounced and the D* tend to small values. When $ increases, so does the importance 
of hydrodynamic interactions between particles and the increasing complexity of particle 
trajectories is reflected in an increasing diffusion constant. When $ increases further, 
then the importance of close particle encounters also increases, restricting the motion and 
consequently reducing the diffusion constants. 

The degree of anisotropy of the diffusion is measured by the ratio Dy/D^, which is 
equal to one in the isotropic case. In the experiments of Nicolai et al. |Q this ratio takes 
values around 5 ... 6 up to $ ~ 0.20. For higher solid fractions the ratio decreases. Such 
a behavior may be expected when the importance of particle encounters increases which 
tends to drive the diffusion towards isotropy. As shown in Fig. ^ we observe a ratio 
Dy/ Dx of about 4 ... 5 at small $ decreasing to ~ 2 ... 3 at large $. 



III.C Hindered settling function 

One of the most easily accessible quantities in a batch settling experiment is the mean sed- 
imentation velocity {V). It may be obtained by observation of the stable upper sedimen- 
tation frontQ, i.e., the interface between batch suspension and clear liquid, or as average 
of the velocity of the bulk particles, if the particles trajectories are traced individually as 
during an index matched settling experiment. The mean sedimentation velocity decreases 
as the solid volume fraction $ increases. For $ — >■ one obtains iy) = Vs. 

The ratio of sedimentation velocity (V) to Stokes velocity is a dimensionless quantity, 
known as the hindered settling function /($), 

/($) = {V)/Vs. (18) 

In [ p6[ an initially step-shaped concentration profile is reported to broaden due to hydrodynamic 
diffusion until a stationary smooth profile is reached after t* ~ 400. 
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The hindered setthng function is always smaller than one, decreases with solid volume 
fraction and is in most low Reynolds number [Re = piVsr /// -C 1] experiments in 3D well 
described by the phenomenological Richardson-Zaki law 0, 

/($) = (!-$)", (19) 

where n typically takes values around 5 in experiments. Theoretical considerations by 
Batchelor et al. lead to 

/($) = ! -5.6$ (20) 

for $ — >^ and Pe ^ 1. Mills et al. p| use mean field and scaling arguments to arrive at 
the expression 

/($) = (l-$)/[l + A;<|./(l -$)='], (21) 

where k = 4.6 in the non-Brownian regime is a parameter chosen to match Batchelor's 
result (|20|) . Thus, these expressions all agree in the low concentration regime and are very 
well established experimentally |2^. Moreover, Eqs. (|T^ and ( pi] ) are accurate also at 



intermediate and high concentrations, not too close to the random packing threshold. 

We display our simulation results in Fig. ^ where /($) is plotted vs. the effective 3D 
volume fraction $ of the simulation. Clearly, {V) decreases when the solid volume fraction 
$ increases. However, the decrease is not as sharp as expected from the Richardson-Zaki 
law (|1^) or from Eq. (|21|). Our 2D simulation results are rather fitted approximately by 
/($) ^ (1 — <l>)'^ for small $. This function is also plotted in Fig. ^ for comparison. 

Apart from the fact that our simulations were done in 2D in contrast to the 3D 
experiments, two factors should lead to a sedimentation velocity that is too high compared 
to the experiments, (i) There is no provision in the simulation to take the volume of 
liquid into account that is displaced by the particles. Thus, our simulation does not 
fully reproduce the net upward backflow effect that otherwise counteracts the settling. 
However, some backflow effect is retained, because the frictional coupling between liquid 
and particles effectively attaches some liquid to every particle. When the particles settle, 
this liquid replaces the fluid present in the lower part of the container and a backflow 
results Q. (ii) The lattice spacing sets the scale on which we resolve details of the 
flow. Particles act on the liquid by point-like drag forces located on the four grid points 
closest to the particle center. The coarser the grid, the less accurate is this model for the 
momentum transfer onto the liquid phase and we underestimate frictional drag because 
we underestimate the interstitial liquid velocities. 

In fact, reducing the grid spacing leads to a more realistic representation of the in- 
terstitial particle velocities and as a consequence to smaller sedimentation velocities. We 
address this question more systematically in Fig. |b, where we plot the simulated values of 
/($), for several values of $ = 0.079, 0.16, 0.24, 0.31 0.39 vs. the ratio Ax/f. In accord 
with the argumentation above, the simulated sedimentation velocity slightly decreases for 
most values as the grid becomes finer and resolves more details of the interstitial flow 
field. 



III.D Velocity fluctuations 

Another interesting quantity is the distribution of velocities of the sedimenting particles. 
In particular, the standard deviation Avy of the velocity distribution has been measured 
by Xue et al. and Nicolai et al. ||^ in different experiments. Here we first focus on the 
dependence of the velocity fluctuations on volume fraction at fixed polydispersity p = 0.04, 
corresponding roughly to the particle size dispersion in the experiment by Nicolai et al. 



P 



To calculate the velocity variance Afj^, in the simulation, we calculate in intervals 
of about At* = 1 the standard deviation of the instantaneous y-velocity distribution of 
selected particles. The selection of particles occurs according to to the same rules as 
used for the determination of diffusion constants (cf. Sec. [1II.B|) . Then we time average 
the standard deviations for times t* > 20, which we found sufficient to avoid significant 
infiuence from the initial transient behavior. 

We then normalize Avy by division with Vs to obtain the dimensionless quantity Av*. 
Figure |^a shows Av* as a function of solid fraction $. For comparison, we have added 
the measured values of Refs. ^ into the plot. The two experimental data sets agree 
at large solid fractions, where they show a pronounced decrease of the fiuctuations. At 
small solid fractions, the data from passes through a maximum, whereas the data of 



25| starts at so high $ that no such maximum is visible. Just as the experimental data. 



also the simulation data shows a decreasing branch for large solid fractions and displays 
a maximum of the fluctuations at intermediate solid fractions, as visible in one of the 
experimental data sets. 

We explain the existence of a maximum in the velocity fluctuations by similar factors 
as have been invoked in the case of the hydrodynamic diffusion constants, which also pass 
through a maximum as functions of $. For increasing volume fraction the fluctuations flrst 
increase due to the more important role of hydrodynamic interactions. For larger volume 
fractions, however, the fluctuations are reduced due to the increasing role of particle 
encounters, which are dissipative in nature and tend to reduce fluctuations. 

In Fig. we plot the ratio Av*/Av*. For large solid fractions our results match 
well the experimental data of Nicolai et al., but do not display the expected pronounced 
maximum for lower volume fractions which is visible in the experimental data. This could 
be due to the small system size limiting the fluctuations and thereby flattening the 
Av* curve at intermediate solid fractions. 



IV Conclusions and Outlook 

We have presented a comparatively simple and fast simulation method to calculate prop- 
erties of particle suspensions. Although our simulations have been performed in 2D and 
involve a set of assumptions that may be considered to be rather crude, namely the ex- 
clusive consideration of translational degrees of freedom and a simplifled local frictional 
coupling of the particulate and liquid phase, we have captured the essential features of 
3D sedimentation experiments, like the solid fraction dependence of (i) the particle self- 
diffusion constants, (ii) the velocity fluctuations and (iii) the hindered settling function. 
Two advantages of our method are that it is (i) not a priori limited to 2D or (ii) to small 
Reynolds numbers. Moreover, more sophistication may enter the drag expression in order 
to improve the quantitative agreement with experiments once the program has been ex- 
tended to 3D. We are currently working on the extension of the method both to 3D and 
to Reynolds numbers around 100. 
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Figure Captions 



Figure |T|: Sedimentation of polydisperse spheres at $ = 0.157: (a) particle trajecto- 
ries in the laboratory frame, (b) particle trajectories in the comoving frame of the 
sedimenting cloud. The particle coordinates have been shifted to start at a com- 
mon origin and the coordinate change due to the mean settling velocity has been 
subtracted. 

Figure |2|: Dimensionless x (Q) aiid y (*) mean square displacements [{xl) = (a;^)/f^ 
and (yl) = (y^)/^^] as functions of dimensionless time t* = tVs/f averaged over 
four runs. The two solid lines indicate a linear approximation to the data in an 
intermediate time regime. The slope of these lines equals twice the dimensionless 
diffusion constants D* and D*. 

Figure |^: (a) Dimensionless particle self- diffusion constants D* (^) and D* (Q) as a 
function of solid fraction $. We have displayed errorbars on selected data points 
to indicate deviations of the mean obtained from n averages over different initial 
configurations [$ = 0.16: n = 6, $ = 0.31: n = 5, all other concentrations: n = 3]. 
(b) Dimensionless ratio Dy/D^ of self-diffusivity parallel (Dy) and perpendicular 
[Dx) to settling direction (Q) and the experimental data from Nicolai et al. P 

(+)■ 

Figure ^: (a) Hindered settling function /($). The graph corresponds to data with a 
fixed ratio of grid spacing to particle size Ax/f = 2.5 (Q) and 1.25 {-k) respectively. 
The solid curve is given by (1 — $)'^. (b) To demonstrate the influence of grid spacing 
on the accuracy of our results, we display the hindered settling function for different 
solid fraction $ as a function of the ratio of grid spacing to average particle radius 
Ax/f. From top to bottom the curves correspond to $ = 0.079, 0.16, 0.24, 0.31 
and 0.39 respectively. 

Figure |5|: (a) Dependence of the width of the velocity distribution Av* on the volume 
fraction. We show simulation data for a polydispersity of p = 0.04 (O)? measured 
data by Nicolai et al. for p ^ 0.04 (+) and Xue et al. |^ for p = 0.01 (x). In all 



cases, the fluctuations Av* and Av* have been made non-dimensional by division 
through the Stokes velocity. Errorbars on selected data points indicate the deviation 
of the mean for averages over different initial conditions, their number being equal 
to the numbers in Fig. |a. (b) Ratio of Av* to Av*. We show simulated data at 
p = 0.04 (O) and the experimental results of Nicolai et al. 
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Table 1: List of fixed simulation parameters. Unless otherwise specified, simulations runs 
have been performed using the parameters in this table. 
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